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Abstract. The diffusion of monovacancies in gold has been studied by 
computer simulation. Multiple jumps have been found to play a central role 
in the atomic dynamics at high temperature, and have been shown to be 
responsible for an upward curvature in the Arrhenius plot of the diffusion co- 
efficient. Appropriate saddle points on the potential energy surface have been 
found, supporting the interpretation of vacancy multiple jumps as distinct 
migration mechanisms. 
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Several face-centered cubic metals, such as Ag, Cu, Au, Ni, and Pt, ex- 
hibit an upward curvature in the Arrhenius plot of the atomic self-diffusion 
coefficient close to the melting point |Q. To explain this effect, three mecha- 
nisms have been proposed @: (a) different types of defects (such as monova- 
cancies and divacancies) contributing to the diffusion process with different 
activation energies 0; (b) an intrinsic temperature dependence of enthalpy 
and entropy changes associated with the formation and migration of mono- 
vacancies Q; (c) contributions from double jumps of monovacancies ||. In 
this Letter we report molecular dynamics simulations in gold aimed at char- 
acterizing qualitatively and quantitatively the effect of correlated multiple 
jumps on the temperature dependence of the diffusion coefficient. Further- 
more, from detailed studies of the potential energy surface, we provide a new 
microscopic interpretation of — at least — the double jump mechanism Q 

It has already been known for some time that correlated jumps occur in 

the diffusion of vacancies in solids close to melting |5|, |6|, |7[], but no statistically 

relevant measurements of this effect have been reported yet. Our results, 

based on the analysis of a large number of events, show that multiple jumps 

are present and very effective in bending the Arrhenius plot near the melting 

point. An interpretation of their microscopic nature has also been needed, 
1 It should be pointed out that our simulations do not address mechanisms (a) and (b) 
above, which therefore remain as further possible factors contributing to the curvature as 
suggested by some experiments. 
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in order to integrate them into the framework of those recent developments 
in the theory of diffusion that explicitly take into account dynamical effects 

Empirical or semi-empirical many-body Hamiltonians developed in the 
last years have proved to be much more realistic for metals than pair poten- 
tials used in the past, in particular to describe surfaces or defects. The system 
we have used in our simulation consists of 255 particles with periodic bound- 
ary conditions, interacting via the "glue" many-body potential 0, which is 
well suited for our study since it reproduces rather well the melting temper- 
ature and the thermal expansion of Au, as well as the T = energetics of a 
monovacancy. We have analyzed all the vacancy jumps in runs performed at 
twelve different temperatures, between 1000K and 1550K, at steps of 50K. 
All the runs are done in the microcanonical ensemble, where temperature 
is defined by time averaging the total kinetic energy. At each temperature, 
the lattice spacing is obtained by an independent zero-pressure simulation of 
a perfect crystal. Since melting occurs near 1350K in this model |9[], runs 
at 1400K, 1450K, 1500K and 1550K refer to an overheated crystal which, 
in the absence of free surfaces, remains stable within the simulation time 
scale (mechanical instability was reached around 1600K). We find these runs 
to be very useful to observe the trends of the diffusion mechanisms at high 
temperatures, even if they cannot be directly compared with experimental 
results. With the exception of the runs at 1000K, 1050K and HOOK, where 
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the jump rate is particularly low, the runs are made long enough to contain 
~ 800 jumps in order to obtain sufficient statistics of rare correlated jumps. 
This requires a run length of 57 ns (8 x 10 6 MD steps) at T = 1150K, and 
progressively shorter runs at higher temperatures. 

From the Arrhenius plot shown in fig. [I] (filled circles), it can already be 
seen that migration of monovacancies alone — the only type of defect present 
in our simulation — produce a measurable curvature. The curve bends just 
below T m and continues with a higher slope into the superheated region. 
The quantity plotted is the atomic migration contribution to the diffusion 
coefficient, defined as 

(u 2 (t)) - « lk ) 

6t ' 1 ' 

where t is the run length, u 2 {t) = \r(t) — r(0) | 2 , ul ulk is the asymptotic value 
of the same quantity for a perfect bulk (not diffusing), independently cal- 
culated and subtracted in order to retain only the diffusive part, and the 
angular brackets indicate average over all the atoms. The atomic diffusion 
coefficient, as usually defined, is D multiplied by the ratio of the concentra- 
tion of vacancies in the system at the given temperature to the concentration 
in our model system (fixed to be 1:256 at all temperatures). 

To determine the presence and role of vacancy multiple jumps, we have 
analyzed all the jump events with respect to both their relative times of 
occurrence and of the atomic displacements involved. In a vacancy multiple 
jump two or more atoms move simultaneously. Each atom moves by only one 
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lattice spacing, d, but the vacancy, as a result of this concerted motion, is 
displaced by one or more lattice spacings, depending on the number of atoms 
involved and their relative directions. In a double jump, for example, if the 
two atoms move in the same direction the vacancy is displaced by 2d , but 
if the two atoms move at an angle of 120° the vacancy will be displaced by 
only d. 

We have first analyzed the jump events in terms of their relative times of 
occurrence. At each time step, the occupation number of each Wigner-Seitz 
cell in the fee lattice is computed. Every change in the occupation numbers 
("event") is recorded, along with the time step at which it occurs. A large 
fraction of the jumps has been found to be followed almost immediately [i.e., 
within a Debye period, corresponding to about 40 MD steps (0.29 ps)] by 
a return jump into the original position. It is very likely that these event 
pairs are not real jumps, but rather occasional excursions of a neighbouring 
atom into the cell containing the vacancy. Such jumps are considered as 
miscounts and eliminated from the analysis. If we consider all the jumps 
followed by a return jump within a time r and calculate the average delay 
5, we observe that when r becomes larger than about 50 steps (0.36 ps), 
5 tends to stabilize around 25 steps, nearly independently of temperature. 
We have therefore taken tr = 25 steps (0.18 ps) as the average delay of 
immediate returns, and removed all the jumps followed by a return within 
3tr = 75 steps. Assuming an exponential distribution, this should eliminate 
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about 95% of these miscounts. After this filtering, the distribution of delays 
between successive jumps at T = 1450K is shown in fig. |2|. For large times 
the semi-logarithmic plot is linear, as expected for uncorrelated jumps. For 
delays of less than 200 steps (1.43 ps), however, strong deviations from lin- 
earity indicate the presence of correlated jumps. The same feature, even if 
not so pronounced, is found also at the lower temperatures, down to 1250K. 
Qualitatively similar results have been found in a simulation study of va- 
cancies in Al ||. Thanks to the much larger statistics of events, we can 
here provide statistically valid information about their contribution to the 
Arrhenius plot, and we also provide an interpretation of their microscopic 
mechanism. 

We have classified all jumps into single, double, triple, etc., by group- 
ing together those events separated in time by less than t m = 200 steps. 
The temperature dependence of the frequencies for the different groups is re- 
ported in fig. H|, referring to the range from 1150K to 1450K [computational 
limitations prevented us to obtain sufficient statistics of multiple jumps for 
T < 1150K, while at 1500K and 1550K (far into the superheated region) 
correlated jumps become so complex to make this analysis impractical to 
do]. The apparent activation energies are not simple multiples of the energy 
for single jumps, as would be expected if they were extracted from a distri- 
bution of uncorrelated jumps. The contribution D$ of single jumps to the 
diffusion coefficient, represented by the open circles in the Arrhenius plot of 
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fig. [I], even exhibits a small downward curvature at the highest temperatures 
(superheated region), in contrast with the upward curvature of the total D. 
From this second step of our analysis, it is clear that the upward curvature 
observed in our simulation is due to correlated jumps. 

We have proceeded in analyzing these correlated events in terms of the an- 
gles between the directions of successive vacancy jumps. Reasonable statistics 
was available only for double jumps. In the fee lattice there are five possible 
angles between the directions of two vectors joining three consecutive lattice 
sites: 0° (collinear), 60°, 90°, 120° and 180° (return jumps). Given a vacancy 
in a specific lattice site, there are m equivalent ways it can perform a given 
double jump, where m = 12, 48, 24, 48, and 12 for 0°, 60°, 90°, 120° and 180° 
double jumps, respectively. For a series of random jumps, the frequency, T, 
divided by the multiplicity, m, should then be uniformly distributed among 
the five angles. This is what we find, within the statistical uncertainty, for 
pairs of jumps with delays r > tm- But for r < tm the distribution is rather 
different, as shown in fig. [|. We can observe some clear trends: a small 
increase above the random ratio for 0°, a definite increase for 60°, a small 
decrease for 90° and a definite decrease for 120° double jumps. 180° double 
jumps also increase, but behave quite erratically with temperature. These 
last events correspond to return jumps of the same atom and do not con- 
tribute to the diffusion. In agreement with the previously given definition, in 
the following we will call double jumps only those events involving two atoms. 
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The apparent activation energies for the separate groups are: 1.9, 1.6, 1.9 and 
1.6 ± 0.2 eV for 0°, 60°, 90° and 120° double jumps, respectively. 

The question now arises if these correlated jumps correspond to dynamical 
effects or to different jump mechanisms 0. The first situation would imply 
the fast successive crossing of two single jump barriers in series, the latter 
the crossing of a different, higher barrier, specific for each double jump. To 
shed some light on this point, we have thoroughly investigated the potential 
energy surface in 3iV-dimensional space |]10"| |, around possible configurations 
for double jump saddle-points. We proceeded by trial and error, first looking 
for watersheds between equilibrium configurations with minimization tech- 
niques, and then looking for extrema in the intermediate regions, always 
letting all the 3N atomic coordinates relax in this process. The calculation 
has been done at the equilibrium density for T = 1350K. Relatively low en- 
ergy paths for double jumps can in fact be found. The corresponding saddles 
are not as simple as the traditional saddle point for single jumps (fig. ||), 
which is quadratic and has only one unstable normal mode, but are instead 
combinations of minima, saddle points and crossroads (saddle points with 
two unstable normal modes), as shown, e.g. for the case of 0° double jumps 
in fig. p. However, the difference in energy between these local structures is 
never larger than 0.02 eV. The energies of the lowest points in the barriers 
for the various double jumps with respect to equilibrium are: 1.93, 1.92, 2.29 
and 2.66 eV for 0°, 60°, 90° and 120° double jumps, respectively. These val- 
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ues are to be compared with the apparent activation energies evaluated from 
MD and reported above. Taking into account the errors in the MD estimates, 
and the fact that such a comparison is valid only at low T, the energies for 
0°, 60° and 90° jumps can be considered rather close. This agreement is an 
indication that double jumps constitute distinct jump mechanisms, corre- 
sponding to the passage over characteristic energy barriers. In fact, we have 
also verified that the distance between the two jumping atoms remains close 
to the first-neighbour distance throughout the jump, in most 0°, 60° and 
90° jumps. The MD estimate of the activation energy for 120° jumps is, on 
the other hand, sensibly lower than the height of the barrier found for this 
mechanism. However, the frequency of 120° jumps is much smaller than the 
other types of double jumps, and these events cannot really be distinguished 
from those that would be expected for a distribution of pairs of successive 
uncorrelated single jumps. 

Finally we note that the anomalous number of return jumps could be 
due to the complicated structure of the saddles, where trajectories can be 
momentarily trapped and then turn back. As shown by fig. [| successful com- 
pletion of a double jump also requires the unstable normal mode to change 
near the saddle point. Moreover, we have found a second family of saddle 
points for single jumps on the sides of the traditional one (fig. |5[). In these, 
the jumping atom is still placed halfway between the two half- vacancies, but, 
instead of lying in the (100) plane, it sits either below or above it by about 
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30% of the nearest-neighbour distance. These lateral saddles have the same 
complicated structure as the saddles for double jumps (fig. ||), with a saddle 
point energy of 1.18 eV, and a local minimum of 1.15 eV, lower than the 
energy of the traditional saddle point for single jumps, 1.34 eV. 

In summary, the migration of monovacancies in gold, as simulated by 
molecular dynamics, exhibits an Arrhenius plot with a pronounced upward 
curvature at high temperature. This curvature is due to multiple jumps, in 
which two or more atoms jump simultaneously. An analysis of the potential 
energy surface reveals the presence of distinct saddle points appropriate for 
double jumps, and supports the interpretation of correlated jumps as distinct 
migration mechanisms. 

Useful comments by Gert Ehrlich, Giulia Galli, Gianni Jacucci, Enrico 
Smargiassi and by a referee are gratefully acknowledged. This work has been 
performed partly under the Progetto Finalizzato CNR "Sistemi informatici 
e calcolo parallelo" and under the U.S. Department of Energy Grant No. 
DEFG02-91ER45439. 
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Figure Captions 

Figure 1. Atomic migration contribution to the tracer diffusion coefficient D 
via monovacancies. Filled circles indicate D as measured directly from 
Eq. (M; diamonds indicate D a \\ as evaluated from the jump frequency 
T, with all the jumps treated as if they were single, uncorrelated jumps; 
and open circles indicate the contribution Ds from only single jumps, d 
is the nearest-neighbour distance and fs = 0.78146 is the geometrical 
correlation factor for single jumps. Apparent activation energies are 
shown. 

Figure 2. Time delays between successive jumps of the vacancy. The bin 
width is 200 steps. 

Figure 3. Frequencies of single (S), double (D), triple (T) and quadruple (Q) 
jumps. Error bars correspond to a standard deviation. The lines, whose 
slopes are given as activation energies, are obtained by a weighted fit. 

Figure 4. Distribution of angles between the jump directions of pairs of suc- 
cessive jumps. The jump frequency T is divided by the geometrical 
multiplicity m (see text). 

Figure 5. Saddle on the 3iV-dimensional potential energy surface for single 
jumps. The plot shows a section of the hypersurface along the directions 
of the unstable normal mode (the jumping mode), rj, and the softest 
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stable normal mode, £. Lengths along r\ and C are in A and the energy 
is reported relatively to the energy of the equilibrium configuration. 
The atomic positions in the saddle point configuration and the atomic 
displacements corresponding to r), projected onto the (100) plane, are 
also shown. Note that r\ and £ are both 3 N- dimensional vectors, but 
their components along the cartesian coordinates of the nonjumping 
atoms are so small that cannot be detected in the figure. 

Figure 6. Saddle for 0° double jumps. Conventions are as in fig. |5|. The 
central configuration is a local minimum. £ is the direction of the softest 
mode. Along the second axis, in order to clearly show the relevant 
structures of the energy surface, we have chosen two different half- 
directions i]_ and rj + for the negative and positive part. This is why 
there is a sharp edge where the two sections meet. The two sketches 
on the left show the atomic displacements in the crossroad positions. 
Direction j]_ corresponds to the first half of the jump: the atom in front 
is slower than the one behind it, so that their separation decreases (as 
also seen by the overlap of the two jumping atoms). Direction r] + 
corresponds to the second half of the jump: the atom in front is faster. 
Completion of a double jump therefore requires a change of the unstable 
normal mode to occur near the central position. 
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